Scaling in a general class of critical random Boolean networks 
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We derive analytically the scaling behavior in the thermodynamic limit of the number of nonfrozen 
and relevant nodes in the most general class of critical Kauffman networks for any number of inputs 
per node, and for any choice of the probability distribution for the Boolean functions. By defining 
and analyzing a stochastic process that determines the frozen core we can prove that the mean 
number of nonfrozen nodes in any critical network with more than one input per node scales with 
the network size A'^ as N'^^^ , with only A''^'^^ nonfrozen nodes having two nonfrozen inputs and the 
number of nonfrozen nodes with more than two inputs being finite in the thermodynamic limit. 
Using these results we can conclude that the mean number of relevant nodes increases for large 
as N'''^^, with only a finite number of relevant nodes having two relevant inputs, and a vanishing 
fraction of nodes having more than three of them. It follows that all relevant components apart 
from a finite number are simple loops, and that the mean number and length of attractors increases 
faster than any power law with network size. 
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I. INTRODUCTION 

Random Boolean networks are often used as generic 
models for the dynamics of complex systems of interact- 
ing entities, such as social and economic networks, neural 
networks, and gene or protein interaction networks 
The simplest and most widely studied of these models 
was introduced in 1969 by Kauffman 2] as a model for 
gene regulation. The system consists of N nodes, each 
of which receives input from K randomly chosen other 
nodes. The network is updated synchronously, the state 
of a node at time step t being a Boolean function of the 
states of the K input nodes at the previous time step, 
t — 1. The Boolean updating functions are randomly as- 
signed to every node in the network, and together with 
the connectivity pattern they define the realization of the 
network. For any initial condition, the network eventu- 
ally settles on a periodic attractor. 

Of special interest are critical networks, which lie at 
the boundary between a frozen phase and a chaotic phase 
[1, In the frozen phase, a perturbation at one node 
propagates during one time step on an average to less 
than one node, and the attractor lengths remain finite in 
the limit N ^ oo. In the chaotic phase, the difference be- 
tween two almost identical states increases exponentially 
fast, because a perturbation propagates on an average to 
more than one node during one time step Whether 
a network is frozen, chaotic, or critical, depends on the 
connectivity K as well as on the weights of the different 
Boolean functions. If these weights are chosen appropri- 
ately, critical networks can be created for any value of 
K. 

The nodes of a critical network can be classified ac- 
cording to their dynamics on an attractor. First, there 
are nodes that are frozen on the same value on every at- 
tractor. Such nodes give a constant input to other nodes 
and are otherwise irrelevant. They form the frozen core 
of the network. Second, there are nodes whose outputs 



go only to irrelevant nodes. Though they may fluctuate, 
they are also classified as irrelevant since they act only 
as slaves to the nodes determining the attractor period. 
Third, the relevant nodes are the nodes whose state is 
not constant and that control at least one relevant node. 
These nodes determine completely the number and pe- 
riod of attractors. If only these nodes and the links be- 
tween them are considered, they form loops with possibly 
additional links and chains of relevant nodes within and 
between loops. The recognition of the relevant elements 
as the only elements influencing the asymptotic dynam- 
ics was an important step in understanding the attractors 
of Kauffman networks. The behavior of the frozen core 
was first studied by Flyvbjerg Q. Then, in an analytical 
study of K — 1 networks Flyvbjerg and Kjaer [7] intro- 
duced the concept of relevant elements. This concept was 
generalized to general critical networks by BastoUa and 
Parisi ^, |9|] ■ They gained insight into the properties of 
the attractors of the critical networks by using numerical 
experiments based on the modular structure of the rele- 
vant elements. Finally, Socolar and Kauffman [l^l found 
numerically that for critical K = 2 networks the mean 
number of nonfrozen nodes scales as Nnj ^ iV^/'^, and 
the mean number of relevant nodes scales as Nrei ~ N^^^. 
The same result is hidden in the analytical work on at- 
tractor numbers by Samuelsson and Troein llj, as was 
shown in [l^] ■ An explicit analytical derivation of these 
and other scaling laws was given in For K = 1, 

these power laws are Nnf N and Nrei ~ N^^'^, since 
there is no frozen core in a -ftT = 1 critical network. 

In this work, we will derive the scaling behavior of 
the number of nonfrozen and of relevant nodes in critical 
Kauffman networks with K > 3. Since the scaling be- 
havior is different for K = 1 and K — 2, one could expect 
that the exponents are generally if-dependent. However, 
we will show that the exponents 2/3 and 1/3 found for 
K = 2 are valid also for larger K and for all possible prob- 
ability distributions of the Boolean functions, as long as 
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the network is critical. We also obtain results for the 
number of nonfrozen nodes with two and more nonfrozen 
inputs, and for the number of relevant nodes with two 
and more relevant inputs. 

The outline of this paper is the following. In the next 
section, we introduce a stochastic process that yields the 
frozen core in K = S networks. The mean-field-theory for 
this process is presented in Section [III[ and an improved 
treatment including fluctuations is presented in Section 
IIV[ yielding the scaling behavior of the number of non- 
frozen nodes in critical networks. The next three sections 
are devoted to special points in parameter space, where 
the stochastic process does not generate all of the frozen 
core. In Sections |V] and IVII those points are considered, 
where the stochastic process gives a smaller frozen core, 
and it is shown that "self- freezing loops" generate the rest 
of the frozen core. In Section [VIIl we consider points in 
parameter space, where the stochastic process does not 
generate any frozen nodes, and where self- freezing loops 
are responsible for all of the frozen core. Finally, in Sec- 
tions IVIIII and IIXI we evaluate the case K > i and the 
scaling behavior of the relevant nodes and attractor prop- 
erties. SectionlXl discusses the implications of our results. 



II. A STOCHASTIC PROCESS THAT LEADS 
TO THE FROZEN CORE 

From now on, we set K — 3 and derive explicitly the 
scaling behavior of the nonfrozen nodes. The generaliza- 
tion to larger K and the scaling behavior of the relevant 
nodes will be discussed later. The first step of the cal- 
culation, which is performed in this section, consist in 
defining a stochastic process that determines the frozen 
core. This process is inspired by the one used in for 
K — 2, however it needed to be modified before it could 
be generalized to larger K. The treatment presented in 
the following is based on the existence of nodes with con- 
stant functions (functions in which the output is fixed 
irrespectively of the input) and it therefore applies to all 
critical models that have a nonzero fraction of constant 
functions. Networks with no constant functions, and in 
particular networks with only canalyzing functions will 
be discussed separately. 

Flyvbjerg Q was the first one to use a dynamical pro- 
cess that starts from the nodes with constant update 
functions and determines iteratively the frozen core. Per- 
forming a mean-field calculation for this process, he could 
identify the critical point. We define in the following a 
process that goes beyond mean-field theory and gives ex- 
act results for the frozen core. We consider the ensemble 
of all networks of size N with a fixed number of nodes 
with constant update functions. All nodes with a con- 
stant update function are certainly part of the frozen 
core. We construct the frozen core by determining step- 
wise all those nodes that become frozen due to the in- 
fluence of a frozen node. In the language of this 
process determines the "clamped" nodes. 



In a = 3 network, each node has 3 inputs, and 
there are consequently 2^ — 256 possible Boolean func- 
tions. In order to specify a model, one has to specify 
the probabilities for a node to choose each of these func- 
tions. Instead of performing the calculation in terms of 
all these parameters, it turns out that three parameters 
are sufficient. For the K = 2 networks, we introduced 

3 parameters corresponding to the occurrence of three 
types of Boolean functions. For larger K, there are more 
types of Boolean functions, and we use therefore a differ- 
ent set of parameters. The first parameter is /?, which is 
the proportion of nonfrozen nodes in the network. 1 — (3 
is therefore the proportion of nodes with a constant up- 
date function. We require /3 < 1 for the calculation per- 
formed in this and the following section. The case (3=1 
will be discussed further below. The second parameter 
is L02, which is the probability that a randomly chosen 
node that does not have a constant update function will 
become a frozen node when one of its 3 inputs is con- 
nected to a frozen node. If one input of a node is fixed at 
some value, the node has effectively two inputs left. We 
now consider those nodes that have not become frozen by 
fixing one input, i.e. we are considering the proportion 
1 — 072 of all nonfrozen nodes. The parameter u>i is then 
the probability that such a node becomes frozen when one 
of the remaining two inputs is connected to a frozen node. 
This probability can again be expressed in terms of the 
probabilities of the different possible update functions. 
Thus all the networks with the same parameters 0^2, wi 
and (3 will be treated as of the same type. As we will 
see below, the properties we are interested in will be the 
same not only for the functions that belong to the same 
type of the network (i.e., that have the same parameters 
but possibly different Boolean functions) but also for the 
different types as long as their parameters are such that 
the network satisfies the criticality condition ([3]) derived 
below. This means that we can have critical networks 
with all possible choices of Boolean functions and that 
they will all be characterized by the same exponents as 
a consequence of being critical. 

Now, let us define the stochastic process that deter- 
mines the frozen core. For this purpose, we differentiate 

4 types of nodes, the numbers of which will change dur- 
ing the process, and we place these nodes in 4 different 
"containers" . Initially, all nodes with constant functions 
are placed in a container labelled T ^ and the remaining 
nodes in a container labelled As. In this container are 
all those nodes, for which we do not yet know if they 
are connected to a frozen node. The other two contain- 
ers, labelled N2 and Ni, are initially empty. They will 
contain nodes with one and two frozen inputs that are 
themselves not (yet) frozen. Since the number of nodes 
in the different containers is going to change during our 
stochastic process, we denote the initial values of num- 
bers of nodes in the containers as iV™*, N^^=Nl™ = 
and iV™*, and the total number of nodes as N™^ (this 
is the actual number of nodes in the network). The con- 
tents of the containers will change with time. The "time" 
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wc arc defining here is not the real time for the dynamics 
of the system. Instead, it is the time scale for a stochastic 
process that we use to determine the frozen core. During 
one time step, we choose one node from the container T 
and determine the influence of this node on the nodes 
connected to it. After determining its influence we will 
remove it from the system, and the number of nodes N 
in the system is reduced by 1. Now, for each nonfrozen 
node in container A3 wc ask whether it receives inpiit 
from the chosen frozen node. If this is the case it freezes 
with probability 0^2 due to the influence of this node and 
moves to container T . With probability 1 — ti)2 it docs 
not become frozen and moves to container M-2.. In one 
time step, we therefore move each node of container Mz 
with probability 2>loiIN to the container and with 
probability 2>{\ — ijj-i) j N to the container Mi. Similarly, a 
node from the container Mi receives input from the cho- 
sen frozen node with probability 2/A'', and it will then 
become frozen with probability Wi and will be placed in 
the container T . If it docs not freeze, we place it in 
container A/i, where we find all those nodes that have 
two inputs from frozen nodes and are not frozen. When 
nodes from this container choose a frozen node as an 
input, they automaticly become frozen. During this pro- 
cess, the probabilities loi and io\ will not change since 
the nodes from containers As and A/2, for which we are 
in every time step determining whether they are going 
to freeze, are chosen at random, and moving them from 
the containers will not change probability distribution 
of the functions of the nodes left in the containers. In 
the next time step, we choose another frozen node from 
container T and determine its effect on the other nodes. 
Some nodes move again to a different container, and the 
chosen frozen node is removed from the system. We re- 
peat this procedure until we can not continue because 
either container T is empty, or because all the other con- 
tainers are empty. If container T becomes empty, we 
are left with the nonfrozen nodes. We shall see below 
that most of the remaining nodes are in container A/i, 
with the proportion of nodes left in containers A/2 and 
As vanishing in the limit A?^*"* — > c». If all containers 
apart from container T are empty at the end, the entire 
network becomes frozen. This means that the dynamics 
of the network goes to the same fixed point for all initial 
conditions. 



III. MEAN FIELD APPROXIMATION AND 
THE CRITICALITY CONDITION 



Let us first describe this process by deterministic equa- 
tions that neglect fluctuations around the average change 
of the number of nodes in the different containers. As 
long as all containers contain large numbers of nodes, 
these fliictuations are negligible, and the deterministic 
description is appropriate. The average change of the 



node numbers in the containers during one time step is 

AiVg = 



37V3 

" N 



(1) 



= -1 
AAT = -1 



■ 0^2- 



N 



The total number of nodes in the containers, N, can be 
used instead of the time variable, since it decreases by 
one during each step. The equation for A^3 can then be 
solved by going from a difference equation to a differential 
equation. 



AN3 ^ dN3 

AN ~ 'dN 
which has the solution 



37V3 
' N 



N3 = N' 



3 ^3™ 



where (3 = ^^jj . Similarly, we find 
^2 = 3(1 -^2)]^^' -3(1 -^2) 



(J\finiy 

m = 3(l-a;i)(l-a;2)/37V-6(l-a;i)(l-a;2);^iV' 



-F3(l- wi)(l -W2)- 



/3 



Nf = {l-3{l-iOi){l-uj2)f3)N 



+3(l-2^i)(l-'-2)^iV^ 
+(3a;r(l-a;2)-l)^^Ar3 



(2) 



When l-3(l-a;i)(l-cj2)/3 < 0, the equation AT/ = 0, 
which represents the stopping condition for the process, 
has a solution for an nonzero value N. This sohition 
shows that the number of nonfrozen nodes in each con- 
tainer is proportional to A?"'"'. This means that on an 
average a nonfrozen node has more than one nonfrozen 
input. A perturbation at one node propagates during one 
time step on an average to more than one node and we 
are obviously in the chaotic phase. 

For 1 - 3(1 - wi)(l - 0J2)P > the equation Nf = 
does not have a nonzero solution for G [0,Af™']. In 
this case, we will stop the process when Nf drops below 1. 
We are in the frozen phase, or we have a critical system. 

In the case 1 — 3(1 — t^i)(l — lji)^ > 0, the values A^3 
a nd A"2 will sink below 1 when A^ becomes of the order 
VN^™, and the higher-order terms contributing to Nf 
and A^i can be neglected compared to the first one. For 
smaller N, only frozen nodes and nodes with one input 
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are left. When Nf falls below 1, there remain only a 
constant number of the nodes of type Mi , 



1 - 3{1 - UJi){l - UJ2)P 



The network is essentially frozen, with only a finite num- 



ber of nonfrozen nodes in the limit N''' 



If we 



now choose the inputs for these nodes, we obtain simple 
loops with trees rooted in the loops. This property of the 
frozen phase was also found in [lO|. 

When parameters of the networks are such that 



l-3(l~c^i)(l-c^2)/3 = 



(3) 



is fulfilled, we are at the boundary between frozen and 
chaotic phase in the parameter space. Thus the network 
is critical. Since the stochastic process stops at Nf = 1, 
we have 



d\2 



1 



(1 - 2uji) (iV""'^) 

(1 -CJl) 



(l-^l) 



f3 



d\3 



In the limit TV"" — > oo the first term is dominant and the 
number of nonfrozen nodes would scale with the square 
root of the network size if the deterministic approxima- 
tion to the stochastic process was exact. We shall see 
below that including fluctuations changes the exponent 
from 1/2 to 2/3. The final number of A/2-nodes for the 
deterministic process for the critical networks is indepen- 
dent of network size, and the final number of A/s-nodes 
is ~ (iV"")~-^/^ and vanishes for A^*"* 00 . We shall 
see below that the fluctuations change these two results 
to ^2 ~ (iV™)i/3 and N3 - const. 

The deterministic description of our process gives the 
wrong scaling of the number of nonfrozen nodes in the 
case of critical networks, but a correct criticality condi- 
tion ([3]). We are interested in the dynamical behavior of 
the networks in the critical phase and we will from now 
on study only networks with the parameters such that 
the criticality condition 1 — 3(1 — uji){l — W2)/3 = is 
fulfilled. 

Before we proceed by introducing the noise into the 
deterministic equations, there is one more piece of infor- 
mation we can extract from the deterministic description 
of the critical process that is going to help us later in de- 
termining the noise term. Introducing n = N/N^^^ and 
Uj = Nj/N^""^ for j = /, 1,2,3, equations ^ simplify to 
(using the criticality condition) 



ni 



1—0^1 



(n^ - n^) 



1 — 2uj-[ 9 

1 — LOl 



UJl 



1 — LOl 



(4) 



This means that our stochastic process remains invari- 
ant (in the deterministic approximation) when the initial 



number of nodes in the containers and the time unit are 
all multiplied by the same factor. For small n, the ma- 
jority of nodes are in container A/i , since ni = n — 0{n^). 
Now, if we choose a sufficiently large TV*"* , n reaches any 
given small value while Nf ~ ■n?N™'^ is still large enough 
for a deterministic description. We can therefore assume 
that for sufficiently large networks Nf/N ~ n becomes 
small before the effect of the noise becomes important. 
This assumption will simplify our calculations below. 



IV. THE EFFECT OF FLUCTUATIONS 

The number of nodes in container A/) , j = 1,2,3, 
that choose a given frozen node as an input is Poisson 
distributed with a mean jNj/N and a variance jNj/N. 
We now assume that n is small at the moment where 
the noise becomes important, i.e., that the variance of 
the three noise terms is Ni/N = ni/n — 1 — 2n + n'^ = 
1 - 0{n) and 2A^2/A^ = 2n2/n = j^in - n^) = 0{n) 
and 3N3/N = 3/3n'^ = 0{in?). All three noise terms 
occur in the equation for Nf, and since the first term 
dominates for small n, we consider only this term in the 
equation for Nf. In the equations for A^i and A^2, the 
noise term is much smaller than the number of nodes in 
these containers and can therefore be dropped. 

The effect of the noise on the final value of N^ can be 
obtained by the following consideration: as we will see be- 
low, the mean final value of A^3 will be a constant, which 
is independent of A^™*. This means that each node that 
is initially in the container As has a probability of the 
order XjN™ of never choosing a frozen input during the 
stochastic process, and this probability is independent for 
each node. From this follows that the final number N^ 
is Poisson distributed with a variance that is identical to 
the mean. This variance is finite in the limit A^*"* 00 
and it does not affect the final value of N2 or A^i . Since 
we have obtained the variance of the final value of N3 by 
this simple argument, we will not explicitly consider the 
noise term in the equation for A^3. 

We therefore obtain the stochastic version of equations 
([T]), where we need to retain only the noise term in the 
equation for A^^: 



AA^3 = 

AA^2 = 

ANf = 

AN = 



1 



N. 



3iV3 
N 
2Af2 

W ^ /3(1 - ^1)77" 

A^i N2 

-IH -+ 2loi — 

N N 

~1 . 



3 - 



1 



N ^ 
(5) 



The random variable ^ has zero mean and unit variance. 
As long as the rij change little during one time step, we 
can summarize a large number T of time steps into one 
effective time step, with the noise becoming Gaussian 
distributed with zero mean and variance T. Exactly the 
same process would result if we summarized T time steps 
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of a process with Gaussian noise of unit variance. For 
this reason, we can choose the random variable ^ to be 
Gaussian distributed with unit variance. 

Compared to the deterministic case, the equations for 
iVa and N2 are unchanged. Inserting the solution for 
and N2 into the equation for Nf, we obtain 



dNf 



N 



N 



(6) 

with the step size dN — 1 and (^^) = 1. (In the con- 
tinuum limit dN — s- the noise correlation becomes 
{i{N)i{N')) =5{N~N')). This is a Langevin-equation, 
and the corresponding Fokker-Planck-equation is 



dP 
dN 



d 



dN 



N 



l-2wi N 



1 



LUl 



1—0^1 

1 d^P 



oji TV*"* 
N 



p 



2dN] 



(7) 



Since we are investigating networks in the thermody- 
namic limit, keeping only the leading terms will give a 
good approximation. Thus, we can neglect the last term 
in the expression under the partial derivative with re- 
spect to Nf once N/N™"^ has become sufficiently small. 
We are left with the Fokker-Planck equation of the same 



type as the one already studied in 
ent coefficient. 



but with a differ- 



dP 
'dN 



d 



dN 



f 



Nf_ 

N 



HN 



^ 2 dNj 



where /i = (1 — 2cl'i)/(1 — uJi) 
We introduce the variables 



N^ 

Vn 



and y — 



N 



(iV™/^)2/3 



(8) 



(9) 



and the function f{x, y) = (iV™/^)^F(iV/, N). The free 
parameter 7 will be fixed below by the condition that the 
probability distribution of the number of nonfrozen nodes 
is normalized. The Fokker-Planck equation then becomes 



^a^ + -^+l2 



3/2\ 9/ ^ 1 _ Q 
^ ) dx 2 dx"^ 



7^ = 0- (10) 



Let W{N) denote the probability that N nodes are left 
at the moment where Nf reaches the value zero. It is 



W{N) = 
Consequently 

W{N) = 



P{Nf,N)dNf- P{Nf,N-l)dNf 



^ I PiNf^N)dNf 



(N' 
{N' 



f{x,y)dx 



-7-1/3 



G{y) 



with a scaling function G{y). W{N) must be a normal- 
ized function. 



W{N)dN = (Ar™/^)-7-i/3+2/3 / G{y)dy = 1 
Jo 



This gives 7 = 1/3. This condition is independent of the 
parameters of the model, and therefore G[y) and /(x, y) 
are independent of them, too. Now, we have 

W{N) = (iV™V^)-2/3G(y) 

The mean number of nonfrozen nodes is 

N= NW{N)dN ^{N'^^'/^if'^ G{y)ydy, 
Jo Jq 

which is proportional to (A^*"*)2/3. 

The probability W2{N2) that N2 nodes are left in con- 
tainer J\f2 at the moment where container J- becomes 
empty, is obtained from the relation 



N2 = 



1 N' 



1 



1-uji TV*"* 1 - (AT'"') 2 ■ 



Since W{N)dN = W2{N2)dN2, we find that the mean 
number of nodes left in container J\f2 is 

/>oo 

N2 = / N2W2{N2)dN2 ^ / N2W{N)dN 
'0 Jo 

^ '-j-ini\l/3 I „,2, 



y^G{y)dy ^ {N"^^f/\ 



M Jo 



In the same manner we find for the number of nodes left 
in container As 

/>oo 1*00 

N3 = N3W3{N3)dN3 = / N3W{N)dN 
Jo Jo 

\2 /-oo 



y G{y)dy ^ const . 



Thus, we have shown that the number of nonfrozen 
nodes scales with network size TV™* as (A^™*)2/3^ with 
most of these nodes receiving only one input from other 
nonfrozen nodes. The number of nonfrozen nodes with 
two nonfrozen inputs scales as (iV*"*)^/'^ and the number 
of nodes with three such inputs is independent of the 
network size. 



V. SPECIAL POINTS AND CANALYZING 
FUNCTIONS 

For uJi — 1/2, the second term in the Langevin Equa- 
tion ([6]) is zero. In this case the next order term has to 
be taken into account since it is the leading one now. We 
will see that the mechanism of creating the frozen core is 
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different for such systems, but in the end we will find the 
same scaling behavior of the number of nonfrozen nodes. 

Now we have to consider the modified Langevin equa- 
tion 

dNf Nf , r.. ( ^ \^ . 

^ = — + 21-/3)— +£ 11 

and the corresponding Fokker- Planck equation 
(f + 2(l-/3)f-^VV + 

V V J\fmi J j 



dP _ _d_ / 
We again introduce new variables 



1 d^P 



2 dNf ' 
(12) 



( i\Tini\2 \ 4/5 

' ^ (13) 



The 



/TV ^ V2(l-/3) 

and the function f{x,y) = (^(^f^y P{Nf,N). 
Fokker-Planck equation then becomes 

For the probability that N nodes are left when Nf 
reaches zero we obtain 

/ (Arini\2 \ -2/5 ^ 

WiN) . (1^) Giy) 

with a new scaling function G. We have used the fact 
that this probability has to be normalized, which gives 
7 = 1/5. 

Using this result, we find for the mean number of non- 
frozen nodes 



1 / (iV™*) 



2 \ 2/5 ~oo 



/3) 



Giy)dy 
(14) 



For the mean number of nonfrozen nodes left in con- 
tainers A/2 and As we find 



N2 = / N2W2{N2)dN2 = / N2W{N)dN 



^^mi^3/5 

(2(1-/3))^/^ 7o 



(2(1-/3))^/^ 7o 



yi/2G(y)dy 
yG{y)dy 



(15) 



and 



iV3W3(^3)rfA^3 = / A^3W^(^)rfAf 



/3 (iV*™)^/^ Z""" 



2 (2(1 - /3))6/5 7o 



yG{y)dy 



(16) 



We see that the number of nodes which become frozen 
due to the influence of the constant functions is smaller 
than in the case of other critical networks. When we look 
at the parameters for these networks more closely, we see 
that these networks are effectively canalyzing with two 
inputs per node. The probability that a node with two 
inputs is going to freeze during one time step is wi = 1/2 
and this means that the network has Boolean functions 
such that nodes with two nonfrozen inputs effectively be- 
long to the Ci or C2 class of Boolean functions with two 
variables, i.e., canalyzing functions. The class Ci contains 
those functions that depend only on one of the two vari- 
ables, but not on the other one. The class C2 contains the 
remaining canalyzing functions, where one state of each 
input fixes the output. It has been shown in [l^| that in 
K = 2 networks with only this type of functions another 
mechanism of creating the frozen core arises. The only 
condition for this is that the number of nodes from class 
C2 is large enough. We will show that it is exactly what 
happens in the networks we are analyzing now. The num- 
ber of nonfrozen nodes with two inputs and canalyzing C2 
functions is here large enough to allow for the creation 
of the self-freezing loops that are going to increase the 
number of frozen nodes and thus change the scaling of 
the nonfrozen nodes from (TV™')-*/^ to (7V™')2/3. 



VI. CREATING SELF-FREEZING LOOPS AND 
THEIR EFFECT 

We are now considering a reduced network consisting 
of those nodes that are not frozen through the influence 
of the nodes with constant functions. The size of this 
network is iV ~ (7V'"')^/^, most of the nodes have one 
nonfrozen input, N2 — (iV™*)"^/^ have two, and N3 ~ 
(7V™*)^/^ have three nonfrozen inputs. Nodes with two 
nonfrozen inputs have a probability to freeze uji = 1/2 
and as such effectively have canalyzing Boolean functions 
of two arguments, belonging to Ci or C2 class. So, the 
number of nodes with two nonfrozen inputs that belong 
to the C2 class has to be ~ (iV™*)'^/^ as it is the fraction 
of all nonfrozen nodes with two inputs. 

Let us now assume that there exist groups of nodes 
that fix each other's value and do not respond to changes 
in nodes outside this group. The simplest example of 
such a group is a loop of C2 nodes where each node cana- 
lyzes (fixes) the state of its successor once it settles on its 
majority bit (the one occurring 3 times in its update func- 
tion table). These loops, introduced in [l3], are called 
self-freezing loops. They can also contain chains of nodes 
with one nonfrozen input or with two nonfrozen inputs 
and a Ci function between C2 nodes. If a chain between 
two C2 nodes as a whole inverts the state of the first C2 
node, the inverted majority bit of the first C2 node has 
to canalyze the second C2 node. The only effect of nodes 
with Ci functions and those with one nonfrozen input in 
such loops is to delay the signal propagation between two 
adjacent C2 nodes. The procedure of finding self-freezing 
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loops is explained in details in . The number of nodes 
on self-freezing loops is there found by mapping the prob- 
lem of finding a self-freezing loop in a C2 network onto the 
problem of finding the relevant nodes sitting on relevant 
loops in a critical network that contains no canalyzing 
functions at all, but only reversible (where the output 
is changed whenever one of the inputs is changed) and 
constant functions. Using results for these reversible net- 
works obtained in [l3| it was found that the number of 
nodes on self freezing loops scales as ~ iV^/^ where N is 
the number of C2 nodes. 

Obviously, nodes depending on or canalyzed by the 
frozen nodes of the self-freezing loops freeze also, and 
such nodes may lead to the freezing of further nodes, etc. 
We can introduce a dynamical process in order to de- 
termine the total number of nodes that become frozen 
because of the self-freezing loops. This process is almost 
the same as the one we have used for identifying the 
influence of the constant functions on the networks dy- 
namics. We again have four containers where the nodes 
left after determining the influence of the nodes with con- 
stant functions are placed. Initially nodes found to be on 
the self-freezing loops are going to be moved from the 
container with nodes with two inputs, J\^2^ to the con- 
tainer J- . Thus the initial number of nodes in the con- 
tainers is going to be iV^ = ((Ar™)3/5)i/3 ^ (^«m)i/5^ 

AtO = (iv™)3/5 - TVj! ~ (iV™i)3/5 and = (Ar™)2/5^ 

and the total number of nodes is = (A^*"')^/^. Now 
we run the same dynamical process as before determin- 
ing influence of the nodes from the frozen loops on the 
rest of this reduced network one by one and then remov- 
ing them from the system. At the end of this process we 
will again have nodes in the container A2. They can now 
make new self-freezing loops made of C2 nodes with the 

chains of nodes with one nonfrozen input between them. 

1/3 

We can then again move nodes that are on the new 
self-freezing loops to the container !F and run the same 
process again. We can even take over the values of Ni , N2 



and iVs and N at the end of the first process, since N; 



1/3 

frozen nodes moved from container J\f2 are negligible in 
comparison to N2. These processes can be repeated as 
long as the number of nodes of type C2 is large enough 
to allow for the creation of self- freezing loops. The equa- 
tions for the change of A^3 and N2 nodes 



ANa = - 
AN2 = - 



3N3 
N 

2N2 
N 



2N3 

13 N 



values of N, N2 and TVs, found in Equations ([14]), (fTSl) 
and pB]) . as initial values of the variables, these differen- 
tial equations have the solution 



N2 



ZLa.3 



Ni + mwi 2 

(iV0)2 



2^3° 



(18) 
(19) 



The number of remaining J\fi nodes increases in the sec- 
ond process, the number of C2 (those in container A2) 
nodes decreases, thus leading to an increasing weight of 
Afi nodes in the nonfrozen network. 

The repeated process of identifying generalized self- 
freezing loops and the nodes frozen by them breaks down 
when the remaining nonfrozen nodes cannot be consid- 
ered as an effective C2 network any more. This happens 
when in the process of creating self-freezing loops the 
probability that a C2 node is going to be attached to 
the end of the chain of nodes with one nonfrozen input 
(thus making closing self-freezing loop possible) becomes 
of the same order of magnitude as the probability that 
this chain becomes a loop. Since the mean size of the 
loops of nodes with one input is found to be of the or- 
der of \/N [15] the assembly of self-freezing loop becomes 
improbable when N2 ~ VlV. 

This condition gives to leading order 



(20) 



(17) 



or ^ (7V*"«)2/3^ -y^g again have the same scaling of 
the number of nonfrozen nodes with the network size. 
The scaling of the number of nonfrozen nodes with two 
and three nonfrozen inputs with the network size we find 
from (HH) and (HH) to be N2 ~ (Af™^)i/3 and N3 ~ const. 
This is the same scaling we have for the case of all other 
critical networks investigated until now. 

When finding the number of nodes on the self-freezing 
loops and defining our second process we assumed that 
there the influence of the nodes with three nonfrozen in- 
puts per node is negligible. We can check if our assump- 
tion was justified. In the beginning of this process the 
number of nodes with three inputs was ~ (7V™*)2/5^ 
The number of nodes that are initially on self-freezing 
loops is (A^2 )^^'^ — (A^™*)^/^. The mean number of nodes 
with three inputs on the self-freezing loops is then 



apply together to all the successive processes of freez- 
ing the network through the influence of nodes of the 
self-freezing loops. Between each two of them the new 
self-freezing loops have been found and moved from the 
container with N2 nodes allowing for the new process to 
start. The equation for N is AA^ = —1, as before. The 
solution of these equations is obtained by going to dif- 
ferential equations for dN2/dN and dN^/dN. Using the 



(A0)V3^ = con.i. 



In the limit of large network size, only a few (if any) 
self-freezing loops are destroyed by nodes with three non- 
frozen inputs, and this does not change the scaling be- 
havior of the number of nodes on self-freezing loops. 
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VII. NETWORKS WITHOUT CONSTANT 
FUNCTIONS 

A. Case uj^ = 1/2, c^2 = 1/3 

Until now, we have assumed that the network has 
nodes with constant functions. In this section, we con- 
sider networks without constant functions, i.e., with P = 
1. The criticality condition ^ then becomes 

3(1-C^l)(l - 1. 

Although the criticality condition was derived under the 
assumption that the network has a nonvanishing propor- 
tion of frozen nodes (i.e., that < 1), it can be extended 
to P = 1, since it is valid for any (3 arbitrarily close to 
1. Furthermore, decreasing (3 slightly for fixed lui and 
UJ2 moves the system to the frozen phase, indicating that 
a system satisfying the criticality condition with (3—1 
is at the boundary of the frozen phase. As we will see, 
the value of the parameters in the critical networks with- 
out constant functions we are considering here is allowing 
the formation of the self-freezing loops and leads to the 
frozen core of the same size as for all the other critical 
networks. Canalyzing networks and threshold networks 
are examples of this category of networks, and they are 
considered important for biological applications. 

The procedure of creating self-freezing loops in the case 
of networks with nodes with two nonfrozen inputs was in- 
troduced and explained in details in [14]. It is the same 
procedure we have used in the previous section. Using a 
similar line of arguments we can explain the assembly of 
the self-freezing loops for the networks with three inputs 
per node determined with parameters being toi = 1/2, 
LU2 = 1/3 and (j — 1. In this case there is a mapping of 
the problem of finding the nodes on the self-freezing loops 
in this network onto the problem of finding the relevant 
nodes on relevant loops in critical network with three in- 
puts per node and only reversible and constant functions, 
i.e., with uji = UJ2 — and /? — 1/3. Self-freezing loops 
are found by starting with a node and keeping track of the 
connection to those inputs that are able to canalyze this 
node if they are canalyzed themselves. This procedure is 
iterated for these input nodes etc., until a loop is formed 
or until it has to stop because no canalyzing inputs are 
found. Similarly, relevant loops in a critical network with 
LUi — u;2 = are found by starting with a node and keep- 
ing track of the connection to those inputs that do not 
have a constant function. This procedure is iterated for 
the nonfrozen inputs etc., until a loop is formed or until 
it has to stop because no nonfrozen inputs are found. In 
both cases, a connection to an input is made with proba- 
bility 1/3, showing that the two processes can be mapped 
on each other. As we will show in section 9 below, in 
critical networks with three inputs per node and nonzero 
fraction of frozen nodes the number of relevant nodes on 
relevant loops scales as (iV™*)-'^/^. Therefore, we con- 
clude that in the network with = 0^2 = 0, the number 
of nodes on self- freezing loops scales also as (iV™*)^/'^. 



We can now proceed just as in the previous section, but 
with (3 — 1 and = TV™*. We continue making self- 
freezing loops and determining which nodes are frozen 
by them, until N2 ^fN . Inserting this condition in 
Equation (fTQ]) . we find to leading order 

iV3/2 

2 = 1 

leading again to ^ (A™*)^/^. 

B. General case 

Now, let us turn to the case /? = 1 with uji < 1/2. (The 
situation wi > 1/2 is not possible for nonfrozen Boolean 
functions with two inputs.) The probability that a node 
we don't know anything about freezes when connected 
to a frozen node is now W2 > 1/3. Every node has three 
inputs and this frozen node could be any of them. This 
means that on an average a node can be frozen by more 
than one input, and the self-freezing components we look 
for in the network here consist of at least as many nodes 
as those in the previous subsection. However, we do not 
need to know the exact number of frozen nodes in these 
components. We will build only one self-freezing loop 
and move its (A™*)^/-^ nodes to the container J^. Then 
we start the calculation of Section [ll] by setting /3 = 1 — 
(A™*)~2/^. Since uji < 1/2, the leading-order terms of 
the calculation performed in Section [TT] are retained in 
this case, and we can take over all the main results of 
that section. In particular, it follows that a single self- 
freezing loop is sufficient to generate the entire frozen 
core, and we do not need to identify other self-freezing 
loops. As before, the number of nonfrozen nodes scales 
as (A™*)2/3. 

VIII. GENERALIZATION TO LARGER K 

The process introduced in Section |TT] can easily be gen- 
eralized to networks with K > 3. We first consider again 
the case (3 < 1. For network with K inputs we define 
a set of parameters /3 and uJi with i G [1,K — 1]. (3 is 
again fraction of the nonfrozen nodes and uji is the prob- 
ability that a nonfrozen node that has K — i inputs from 
frozen nodes freezes when receiving another frozen input 
in our process. These K parameters are going to define 
completely the class of networks we observe in the pro- 
cess. Using the deterministic description of the process 
analogous to the one described in Section Hill we find the 
criticality condition for networks with any K: 

K{1 - LJl)(l - C^2) • • • (1 - UJK-l) = 1 ■ (21) 

Introduction of noise in the process gives the Langevin 
equation 

dNf Nf , , ^ f N y ^ 

rf/ = ^+E/^(-i'---'-'0(^j +^ (22) 
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where the . . . ,uJi) are functions of the parameters 

of the system obtained from the stochastic process They 
satisfy fi{uji, . . . ,uji) = when ujj = l/(j + 1) for aU 
J G [I,*]- We see that in this general Langevin equation 
the leading term in N is the same as in Equation ([6]). 
Therefore we find that in the thermodynamic limit the 
number of nonfrozcn nodes scales in critical networks as 
(iV*"*)^/^ with the network size. 

Just like in the A' = 3 networks, parameter values can 
be such that one or more of the leading terms in the 
Langevin equation vanish. These special points in the 
parameter space describe networks where the Boolean 
functions are such that the nodes left nonfrozen after de- 
termining the influence of the frozen nodes in our process 
can additionally generate self-freezing loops. Their influ- 
ence on the rest of the network has to be determined by 
generalization of the process introduced in Section IVII 
The number of classes of special points will increase with 
K, leading to a hierarchy of special points. For each K, 
there are K — 3 classes of points in parameter space that 
are equivalent to the special points of networks with K—1 
inputs per node (that is they have the same leading term 
in the Langevin equation), and one new class of special 
points where only the last term in the Langevin equation 
(1^^ is nonzero. Furthermore, there is the case /3 = 1. As 
an illustration, in the case K = 4 there are two classes 
of special points for /3 < 1. One of them has wi = 1/2. 
In this case, the influence of the frozen nodes will lead 
to (iV™*)'*/^ nonfrozen nodes. Boolean functions of the 
nodes with 2 nonfrozen inputs and the number of them 
left after the first process are such that self-freezing loops 
are made and their influence will again give (A^'"')^/^ as 
the number of nonfrozen nodes in the network. This case 
can obviously be reduced to the K = 3 network. The 
other class of special points is obtained when the param- 
eters of the network are uji — 1/2 and W2 — 1/3. In this 
case, (iV™*)^/"^ nodes will be left nonfrozen after deter- 
mining the influence of the frozen nodes. One can easily 
show that the creation of self-freezing loops is possible 
and that their influence leads to a number of nonfrozen 
nodes that scales as (A^*"*)2/3 with the network size. 

For general values of K, the K — 2 classes of spe- 
cial points with /3 < 1 are given by the condition luj — 
+ 1) for all j G where i takes for every class 
one of the values from the interval [1, K—2]. This means 
that /i = 0, . . . , /i = in the Langevin equation 
and the term . . . , a;,+i)(iV/iV*™)*+i is the lead- 

ing one. The nodes left nonfrozen after determining the 
influence of the nodes with constant functions scale with 
the network size as 

(jY™)(2*+2)/(2i-H3)_ rpj^^ numbers and 
Boolean functions of the nodes with k G [2, « -I- 1] non- 
frozen inputs are such that they allow for the creation of 
the self-freezing loops, and their influence will for each of 
these special points, i.e., for each i € [1,K — 2], reduce 
the number of nonfrozen nodes to (A^™')2/3. 

For networks without constant functions (that is with 
(3 = 1) the frozen core arises only because of the cre- 
ation of self-freezing loops and their effect on the network. 



Just like for all other parameter values, there is straight- 
forward generalization of the analysis performed for this 
type of networks in the case when K ^ 3 in Section IVlII 
In the case when Ui = l/(i -I- 1) for alH G [1, iiT — 1] there 
exists again a mapping of the self-freezing loops on the 
relevant loops of a X critical network with only reversible 
and nonfrozen functions, from which it follows that the 
number of nodes that are initially on self-freezing loops 
scales as (iV*"*)i/3. The process described in Section IVll 
can then be generalized to these networks. For any other 
choice of parameters satisfying the criticality condition 
(PT|) for (3 = 1, self- freezing loops can also be formed, 
and after moving only one of them in the container with 
frozen nodes we will have the same process as for the one 
of the classes of critical networks with (3 < 1 that were al- 
ready studied. Scaling of the number of nonfrozen nodes 
in the critical networks without frozen nodes and any 
fixed number of inputs will be the same as in all other 
critical networks. 

Let us end this section by noting that there is another 
class of special points when the Boolean functions are 
chosen such that each of them responds only to one of 
the K inputs. In this case, the network is effectively a 
K = 1 network, since for each node those K—1 inputs 
to which the node does not respond, can be cut off. In 
the calculations of the previous sections we have always 
assumed that a nonvanishing proportion of functions is 
not of this type. 



IX. RELEVANT NODES AND THE NUMBER 
AND LENGTH OF ATTRACTORS 

Relevant nodes are the nodes whose state is not con- 
stant and that control at least one relevant node. These 
nodes determine completely the number and period of 
attractors. In the previous sections, we have shown that 
the number of nonfrozen nodes scales as (A^™*)2/3 for 
any critical network. We have also seen that among them 
there are only (iV™*)i/3 nodes having two nonfrozen in- 
puts, and that the number of nonfrozen nodes with more 
than two nonfrozen inputs vanishes in the thermody- 
namic limit. The nonfrozen nodes can now be connected 
to a network. This is a reduced network, where all frozen 
nodes have been cut off. In jl^, we defined a stochas- 
tic process for the formation of this reduced network and 
the identification of the relevant nodes for critical K = 2 
networks. The relevant nodes are determined by remov- 
ing iteratively nodes that are not relevant because they 
influence only frozen and irrelevant nodes. The number 
of relevant nodes was found to scale as (A^™*)i/3^ and 
the scaling function characterizing their probability dis- 
tribution depends on the parameters of the model. 

The scaling of the number of nonfrozen nodes as well 
as the scaling of the number of nonfrozen nodes with 
two nonfrozen inputs as function of the network size is 
the same for every critical network, as we have shown 
in this paper. Since the fraction of nodes with more 
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than two nonfrozen inputs is vanishing in the thermo- 
dynamic hmit, the network of nonfrozen nodes, which 
is the starting point for the process of determining the 
relevant nodes, is the same as in the K — 2 case. So, 
we can conclude that the results for the scaling of the 
number of relevant nodes found in [13] for the K = 2 
critical networks are valid for any critical network. The 
number of relevant nodes in critical networks scales as 
(iV™*)i/3 with the network size. Among them are a con- 
stant number of relevant nodes with two relevant inputs 
and a vanishing number of relevant nodes with more than 
two relevant inputs in the limit iV™* ^ oo. If only these 
nodes and the links between them are considered, they 
form loops with possibly additional links and chains of 
relevant nodes within and between loops. 

It follows that all critical networks with K > 1 show 
the same scaling behavior. The only exception is the case 
_ftr = 1, which is different because there is no frozen core. 

As we have shown in [l3l |. we can derive properties 
of attractors from the results for the relevant nodes. In 
particular, we can take over the result of [l3j that all rel- 
evant components apart from a finite number are simple 
loops, and that the mean number and length of attrac- 
tors increases faster than any power law with the network 
size. 



X. CONCLUSIONS 

In this paper, we have considered the limit of large 
network size, and we have found the scaling behavior of 
the number of nonfrozen nodes, of the number of non- 
frozen nodes with more than one nonfrozen input, of the 
number of relevant nodes, and of the number of relevant 
nodes with more than one relevant input in a general class 
of critical random Boolean networks with fixed number 
of inputs per node. The mean values of these quanti- 
ties scale with network size N"^™ as a power law in TV™' , 
with the exponents being the same for any critical net- 
work. No matter what the distribution of the Boolean 
functions is and how many inputs per node the criti- 
cal network has, number of nonfrozen nodes scales with 
the network size as (A^™')2/3^ ^jig number of nonfrozen 
nodes with two nonfrozen inputs scales as (iV™*)^/'^, the 
exponent for the number of nonfrozen nodes with three 
nonfrozen inputs is zero, and it is —n/3 for the number of 
nonfrozen nodes with n + 3 nonfrozen inputs. The num- 
ber of relevant nodes scales always as (A^™*)i/3^ with a 
constant number of them having two inputs and a van- 
ishing proportion having more than two. 

It follows that all critical random Boolean networks 
with K > 1 belong to the same class of systems. Chang- 
ing the weights of the different Boolean functions (for 
instance by choosing threshold networks or canalyzing 
networks) or changing the number of inputs per node 
(which might make the model more relevant for biological 
applications) will not change the scaling of the number 
of nonfrozen and relevant nodes with the size of the net- 



work, and it will not change the fact that the number and 
length of attractors increases faster than any power law 
with the network size, as long as the network is critical. 
Using a different method, Samuelsson and Socolar have 
recently also found that the number of nonfrozen nodes 
scales in the same way for all > 1 critical networks 

M- 

From the calculations performed in this paper it can 
be concluded that the results are also valid for networks 
that have nodes with different values of K. If Kmax is the 
largest number of inputs occuring in the network, we can 
set K = Kjnax, and we can view nodes with less inputs as 
nodes with Kmax inputs, but with a function that does 
not depend on all of its inputs. In contrast, our results 
cannot be generalized to networks with a broad distri- 
bution of the number of outputs. The method employed 
in this paper is based on a Poissonian distribution of the 
number of outputs, and is most likely valid also for other 
distributions as long as the second moment of the number 
of outputs is finite. This can for instance be concluded 
from the analogy between the propagation of activity in 
a Boolean network and percolation on a directed graph, 
for which many results are known [l3]- 

The finding that the number and length of attractors 
in critical Boolean networks increases superpolynomially 
with network size is detrimental to the hypothesis that 
these networks are models of gene regulation networks, 
where only a limited number of dynamic pathways should 
exist. However, by considering asynchronous update in- 
stead of parallel update and by requiring that dynam- 
ics should be robust with respect to fluctuations in the 
update sequence, the number of attractors reduces to a 
power law in system size, which is more realistic than the 
superpolynomial growth [3 [3 • The method presented 
in this paper is independent of the updating scheme, and 
the scaling of the number of nonfrozen and relevant nodes 
is therefore same for asynchronous update as for paral- 
lel update. The relevant components are consequently 
also the same. With the insights obtained in the present 
paper, we can immediately apply the results for asyn- 
chronous update in = 2 critical networks to critical 
networks with larger values of and we can conclude 
that the number of attractors in critical networks with 
asynchronous update increases as a power law of the sys- 
tem size. 

Finally, let us consider networks where the connections 
between nodes are not made at random, but that show 
some degree of clustering. Such networks have a finite 
proportion of nodes that have correlated inputs and that 
can therefore become frozen, e.g., because their inputs 
are always in the same state. In contrast, the randomly 
wired networks considered in the present paper have only 
a limited and small number of nodes with correlated in- 
puts even in the thermodynamic limit of infinite network 
size. For small-world networks, which have a high de- 
gree of clustering, our method for determining the frozen 
core is not valid, because it is based on the assumption 
that nodes choose their inputs independently from each 
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other. Small- world networks need therefore a separate Acknowledgements: we thank Viktor Kaufman for use- 
analytical treatment, which has not been done so far. ful discussions. 
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